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ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОМЕРНОЙ ПЛОСКОЙ ЗАДАЧИ СТЕФАНА 


Рассмотрен подход к моделированию и численному решению задачи о распространении фронта фазового 
превращения в мерзлых породах, который весьма актуален в связи с бурением и эксплуатацией скважин в 
условиях вечной мерзлоты. Приведено точное решение задачи в случае полубесконечного физического про- 
странства. С помощью вычислительного эксперимента определена область его применения. 
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Введение. Рассматривается одномерная модель распространения плоского фронта фазового 
превращения в твердой среде. Данная задача имеет множество приложений. В частности, она 
может быть рассмотрена в связи с поддержанием в надлежащем тепловом состоянии приустьевой 
зоны скважин, эксплуатирующихся в условиях вечной мерзлоты [1]. Модель основана на следую- 
щих представлениях. В каждый момент времени физическое пространство разделено на две зоны, 
различающиеся тем, что заполняющее их одно и то же твердое вещество находится в разном аг- 
регатном состоянии и, соответственно, обладает разными теплофизическими свойствами. По- 
следние полагаются однородными по пространству и неизменными во времени для каждой из 
зон. Область фазового превращения вещества из одного состояния в другое имеет малую толщи- 
ну, что позволяет считать её поверхностью. При переходе через эту поверхность свойства среды 
меняются «скачком». 

При указанных допущениях математическая модель задачи в безразмерных переменных и 
величинах включает в себя уравнения распространения тепла для каждой из зон: 
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Здесь 27 — пространственная координата; { — время; Й - продольный размер физической области; 
Т - температура в точке с координатой 2 в момент времени & 1, Т, - масштабные значения тем- 
пературы, первое из которых соответствует ее значению при 2=0, а второе - при 2=1; с, д, р - 


соответственно удельная теплоемкость, коэффициент теплопроводности и плотность вещества; 
индекс ф относится к величинам на границе раздела между зонами; 1-к левой зоне, 2 -к пра- 


вой; 2. - Координата фронта фазового превращения. 


Кроме того, модель включает граничные условия первого рода: 
при у=0: ®=0, (2) 
при у=1: 9=1, (3) 
что соответствует тому, что на одной из границ на протяжении всего периода наблюдения под- 
держивается температура 7, а на другой - Т,. 
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Начальное распределение температуры полагается однородным: 
при т=0, уЕ(0;1]:©9 =1. (4) 
Дополнительно на “межфазной” границе требуется выполнение условий непрерывности 
температуры и теплового баланса, то есть при т > 0, у=у,: 








9(7,-0)=00,+0)=9,; (5) 
29| 9 оу е 
ду У6-0 ду Уь+0 
а Т.-Т 
где р 9% вы ЕР Че й и. - температура фазового перехода; 


де РЕ а 





4, - теплота фазового перехода в расчете на единицу массы вещества в состоянии 1; 


Г - безразмерная скорость перемещения фронта фазового превращения. 

Заметим, что задача (1) - (6) соответствует процессу фазового превращения вещества, 
первоначально находившегося в однородном состоянии с температурой Т, , когда на левой гра- 
нице, начиная с начального момента времени и в дальнейшем, устанавливается температура 7’. 
Подразумевается, что Т, < 7, <1,. 

Из (4) в силу очевидного неравенства 0<®„<1 следует, что начальная координата 
фронта у,(0)=0. 

Для расчета начальной стадии процесса, а также для тестирования численного алгоритма 
применим аналитическое решение задачи для полубесконечного пространства. 

Для малых значений переменной времени т влияние правого граничного условия (3) не- 


существенно для развития процесса. Поэтому заменим его условием на бесконечном расстоянии: 
при у=®: 9 =1. (7) 


Будем искать решение задачи (1),(2),(4)-(7) в виде ®(т;у)= Г (1), где п= 


= 


Уравнение (1) удовлетворяет функция: 


п 2 
о +55, при 0О<п<® 
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(т) = . - 
Си 
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К 
где 7, = А, соответственно, у, = к/т ; Г = ——. Чтобы выполнялись граничные условия и ус- 


2/т 
ловия на границе раздела “фаз” (при 77, = К), константы в (8) должны определяться из следую- 


щих выражений: 
Ь, =0, 
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Как видно, приведенные соотношения зависят от безразмерной координаты границы раздела 
П» = К. Последняя, в свою очередь, определяется из трансцендентного уравнения: 


й, =©.. (9) 

В качестве модельных примем следующие значения входящих в задачу констант: 

А, = 1,65 Вт/(м-К); А= 1,32 Вт/(м-К); с,=1220 Дж/(кг.К); с =1640 Дж/(кг.К); 
р. =1780 кг/м*; р, =1800 кг/м? ; 4, =66740 Дж/кг; й=10м; Ти=+6°С; 72 = -29С; Т, =0°С. 

Решением уравнения (9) для базового набора параметров является К = 0,5003. Соответ- 


ственно, зависимость от времени координаты фронта фазового превращения в задаче с полубес- 
конечным физическим пространством определяется безразмерным уравнением: 


ур = 0,5003-/т . (10) 

Для численного решения задачи (1)-(6) с конечными размерами области предварительно 
перейдем к ее дискретному аналогу. Для этого область решения равномерно разобьем на элемен- 
тарные контрольные объемы [/_,;и], {= 1,2.....п, центры которых примем в качестве узловых 


точек у,, и проинтегрируем уравнение (1) по каждому из них [2]. Заметим, что на каждом вре- 


менном шаге один из контрольных объемов будет содержать поверхность фазового превращения. 
При интегрировании по контрольному объему, в котором происходит фазовый переход, учитыва- 
ется условие теплового баланса (6). В результате для расчета температурного поля на каждом 
шаге по времени получаем неявную систему линейных уравнений относительно температур в уз- 
ловых точках с трехдиагональной матрицей: 


&+1 &+1 К-1 
а,0; = с 7 46“: +Ь,, (11) 
1 1 
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5 5 5 Ах 5 5 Аг 
К-1 А А А А 
У =УОНУЮ Аг; у =у,(т,}; Г =У(т,)} та =т,+Ат; 5 - индекс кон- 
трольного объема, в котором происходит фазовое превращение. 
Как видно, кроме температур система (11) включает неизвестное значение /“” скорости 


фронта фазового превращения. Для решения полученной системы применялся метод сквозной 
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«прогонки» в сочетании с методом «пристрелки» в отношении величины У“^, что позволяет на- 
ходить попутно очередную координату поверхности раздела “фаз” у%*”. 


Результаты проведенных на компьютере численных расчетов для базового набора значе- 
ний параметров представлены ниже в графическом виде (рис.1). Абсцисса точки излома на каж- 


дой из кривых соответствует координате у, фронта фазового перехода, а ордината этой точки 


равна температуре фазового превращения, которая для базовых величин имеет значение 
©, = 0,75 . Последнее из представленных на рис.1 распределений температуры близко к стацио- 


нарному распределению, когда скорость распространения поверхности раздела “фаз” равна нулю. 

Для того чтобы рассчитать положение У “стационарного” фронта достаточно в уравнении (1) 

производные по времени положить равными нулю и воспользоваться условиями теплового балан- 
9, 

©,+Л^(-©,)’ 





са (5), (6). В результате получаем соотношение и = 


Температура © 








0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,5 0,9 1 


Координата у 


Рис.1. Распределение температуры @®© в различные моменты времени т: 


1-т= 0,0002;2-т=0,0777;3-т=0,6557;4-т =11,71 


Для базовых параметров это дает У _ 12, 0,706. Численный расчет, что в некоторой 
Е И 
степени иллюстрирует рис.1, действительно приводит к данному значению. Это служит одним из 
подтверждений того, что разработанный численный алгоритм вполне соответствует описываемо- 
му процессу. 
На рис.2 представлены зависимости от безразмерного времени т безразмерных скорости 


Г фронта фазового превращения и координаты самого фронта у,, полученные в результате 


численного решения задачи (1)-(6). Здесь же в виде набора точек показано определяемое из со- 
отношения (10) изменение с течением времени положения границы раздела между зонами для 


задачи с бесконечной протяженностью физического пространства. Как показывают расчеты, 
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влияние правого граничного условия в задаче (1)-(6) несущественно, пока фронт фазового пере- 
хода у, не достигает значений близких к 0,2 (соответственно, пока переменная времени т не 


больше 0,16). Об этом свидетельствуют, в частности, почти идентичные распределения темпера- 
Тур на каждом шаге по времени в совместной расчетной области для каждой из задач. 


Координата фронта, ух 
Скорость фронта, \ 





0 0,2 0,4 0,6 0,8 1 1,2 1,4 


Время, т 


Рис. 2. Изменение поверхности раздела между зонами с различным агрегатным состоянием 
вещества: 1, 2 — соответственно графики скорости перемещения фронта и координаты фронта; 
3 — график изменения с течением времени положения границы раздела “фаз” 


Эти распределения определяются по формулам (8). Кроме того, при численном решении 
задачи (1)-(6) отслеживались изменения величин у, / Л и 2. ‚ которые в задаче с беско- 


нечной толщиной грунта представляют собой константу А. Оценивать относительную величину 
разницы между значениями у„, полученными для каждой из двух моделей, не имеет смысла, по- 


тому что, во-первых, для малых значений времени малы и значения у’, и, в силу того, что ЭВМ 


работает с ограниченным числом разрядов, при малых т относительная разность может быть 
довольно большой, а, во-вторых, начальная скорость перемещения фронта в силу модельного 
представления о разрывности начального распределения температур должна быть бесконечной, 
чего также невозможно добиться при численном решении с помощью ЭВМ. С учетом указанных 
обстоятельств, для проверки близости решений разумным представляется соотнесение величин 


К, = У / и и К, =2иШ ‚ получаемых в ходе численного решения задачи с конечным про- 
дольным размером физического пространства, со значением А, определяемым из (9). 

Заключение. Расчеты показывают, что на начальном временном промежутке относительные 
разницы между К, и К, атакже К, и К, не превышают 10% от А, причем чем больше т отли- 
чается от нуля, тем меньшей становится эта разница, и лишь тогда, когда начинает существенно 


сказываться условие на правой границе задачи (1)-(6) (соответствующие т намного больше зна- 
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чения, которое соответствует концу начального промежутка), она приобретает явную тенденцию 
к возрастанию. Заметим, что при неограниченном возрастании переменной т кривая 2 (см.рис.2) 


(с) 


выходит на горизонтальную асимптоту у= у, (для базовых параметров у = 0,706). Последнее 


означает, что положение фронта стабилизируется. При этом кривая 3 асимптоты не имеет, то 
есть координата поверхности раздела зон в случае бесконечной расчетной области не имеет 
верхнего предела. 
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